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Abstract 


This paper documents the CFL3D contribution to the AIAA Supersonic Shock Boundary Layer 
Interaction Workshop, held in Orlando, Florida in January 2010. CFL3D is a Reynolds-averaged 
Navier-Stokes code. Four shock boundary layer interaction cases are computed using a one-equation 
turbulence model widely used for other aerodynamic problems of interest. Two of the cases have 
experimental data available at the workshop, and two of the cases do not. The effect of grid, flux 
scheme, and thin-layer approximation are investigated. Comparisons are made to the available ex- 
perimental data. All four cases exhibit strong three-dimensional behavior in and near the interaction 
regions, resulting from influences of the tunnel side- walls. 


1 Introduction 

Shock boundary layer interaction (SBLI) is a complex supersonic fluid dynamic phenomenon of 
relevance in many applications of aerodynamic interest, particularly for supersonic inlet flows. The 
shock can cause boundary layer separation, which is typically an unsteady phenomenon. Fig. 1 
shows a schematic diagram of a SBLI in 2-D. The angle 9 of the shock-generating plate defines the 
strength of the oblique shock wave (the larger 0, the stronger the shock). The oblique shock wave 
extends down to the lower wall, where it interacts with the oncoming turbulent boundary layer. The 
computation of SBLI cases has been a focus of the CFD community for many years [1]. 

In January 2010, a SBLI workshop was held in conjunction with the American Institute of Aero- 
nautics and Astronautics (AIAA) 48th Aerospace Sciences Meeting in Orlando, Florida. The pur- 
poses of the workshop were (1) to assess the state-of-the-art computational methods for predicting 
turbulent SBLI flows, (2) to provide an impartial forum for evaluating the effectiveness of current 
computer codes and modeling techniques using Navier-Stokes solvers, and (3) to identify areas 
needing research and development. Several papers came out subsequently, describing the work- 
shop, its results, and lessons learned [2-5]. The reader is referred to these papers for additional 
details regarding the workshop that are not covered here. 

This paper documents a particular submission to the workshop: that of the structured Reynolds- 
averaged Navier-Stokes (RANS) code CFL3D [6]. The first part of the paper briefly describes the 
numerical method. Then, a short description of the experiments is provided. Results are then given 
for four different cases and conclusions are drawn. 


2 Numerical Method 

CFL3D is a multi-zone compressible Navier-Stokes CFD code in wide use in U.S. industry. It uses 
point-matched, patched, or overset grids, and employs local time step scaling, grid sequencing, and 
multigrid to accelerate convergence to steady state. A time-accurate capability is also available, but 
for all the SBLI cases only steady-state simulations were performed. 

CFL3D is a finite volume method. It uses third-order upwind-biased spatial differencing on the 
convective and pressure terms, and second-order differencing on the viscous terms; it is globally 
second-order spatially accurate. Both flux vector splitting (FVS) [7] and flux difference splitting 
(FDS) [8] are available to obtain fluxes at the cell faces. Unless otherwise noted, for all of the 
computations in this paper, FVS was used due to its greater robustness for strong shock cases. The 
code is advanced in time with an implicit three-factor approximate factorization method. For most 
of the applications in this paper, the thin-layer form of the Navier-Stokes equations was solved, with 
viscous terms active in each of the three coordinate directions in computational space (i.e., viscous 
cross-derivative terms were neglected). Recently, full Navier-Stokes capability has been added to 
CFL3D; this capability was exercised for one of the cases to show its effect. 
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Table 1 . Summary of the Experimental Cases 


Organization 

0 , deg. 

Moo 

To, K 

Po, kPa 

Reo 

status 

IUSTI (UFAST) 

8.0 

2.25 

293 

50.5 

6900 

open 

UM (CCAS) 

7.75 

2.75 

293 

101 

6600 

open 

UM (CCAS) 

10.0 

2.75 

293 

101 

6600 

blind 

UM (CCAS) 

12.0 

2.75 

293 

101 

6600 

blind 


For supersonic cases such as SBLI, a reconstruction limiter is required to suppress spurious 
oscillations near the shock. A smooth limiter tuned to the upwind-biased spatial differencing scheme 
was employed, as described in Krist et al. [6]. 

There are several turbulence models available in CFL3D. All are solved loosely-coupled to the 
Navier-Stokes equations. For the current SBFI applications, only the one-equation Spalart-Allmaras 
(SA) model [9] was employed. This turbulence model is widely used throughout the aerospace com- 
munity for applications that include subsonic, transonic, and supersonic flows. It is also considered 
to be a good model for predicting both attached flows as well as flows with small regions of separa- 
tion, including shock-induced separation. For all applications here, the model was active everywhere 
(fully turbulent flow was assumed). 


3 The Experiments 

There were four sets of experiments used for the workshop. One of these was conducted at the 
Institut Universitaire des Systemes Thermiques Industriels (IUSTI) in Marseille, France [10], in 
support of a European Union project titled UFAST [11]. The other three were conducted at the 
University of Michigan (UM) [12], as part of the U.S. Air Force sponsored Collaborative Center for 
the Aeronautical Sciences (CCAS). Two of the cases were open (experimental data were supplied 
to workshop participants ahead of time) and two were blind (no experimental data were available 
to participants). Some of the details regarding the four experiments are given in Table 1. For the 
IUSTI (UFAST) experiment the shock generator spanned the tunnel, whereas for the UM (CCAS) 
experiments it did not. Additional details can be found in DeBonis et al. [3] or on the workshop’s 
website 1 . Particle Image Velocimetry (PIV) was used to obtain measurements of velocities as well 
as turbulence quantities in the interaction regions of the experiments. DeBonis et al. [3] also has a 
good discussion of the experimental uncertainties. 

It is important to note that these SBFI cases are three-dimensional. This realization was one of 
the major outcomes of the UFAST project [4, 13]. Comer vortices in the tunnel test section affect 
the entire flow field, including the shock and interaction region in the center of the test section away 
from the corners. Therefore, computations must include the three-dimensional geometry, including 
all tunnel walls, in order to properly account for the flow physics. Some of the three-dimensional 
behavior of these cases will be shown in the Results section. 


4 Results 

Table 2 summarizes the computer runs performed using CFF3D. The fine grids for each case, con- 
sisting of between 7 and 9 million cells, were created at NASA Glenn Research Center by N. Geor- 
giadis, M. Vyas, and D. Yoder for use by both the WIND code at NASA Glenn and the CFF3D 
code at NASA Fangley. All grids assumed span wise symmetry; i.e., one span wise boundary was 
modeled as a solid wall, while the other spanwise boundary (the centerline plane) used a symmetry 

1 http://sbhworkshop.engin.umich.edu/html/, cited 8/12/2010. 
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Table 2. Summary of CFL3D Cases Run 


Case 

Grid size (millions of cells) 

flux scheme 

thin vs. full 

IUSTI (UFAST) 9 = 8.0° 

7.07 

FVS 

thin 


0.88 

FVS 

thin 


0.11 

FVS 

thin 

UM (CCAS) 0 = 7.75° 

7.84 

FVS 

thin 


0.98 

FVS 

thin 


7.84 

FDS 

thin 


0.98 

FDS 

thin 


0.98 

FDS 

full 

UM (CCAS) 9 = 10.0° 

8.50 

FVS 

thin 

UM (CCAS) 6 = 12.0° 

8.50 

FVS 

thin 


boundary condition. A limited number of grid convergence studies were performed for the IUSTI 
(UFAST) 0 = 8.0° and UM (CCAS) 0 = 7.75° cases, and the effect of flux scheme and thin-layer 
vs. full Navier-Stokes was briefly looked at for the UM (CCAS) 0 = 7.75° case. The results of these 
studies will be described in their respective sections below. 

4.1 UFAST Case 

An overall view of the UFAST fine grid is shown in Fig. 2. All distance units are in mm. This grid 
contained 10 zones, most of which were connected in a point-matched (one-to-one) fashion. There 
was one patched interface located between the two zones approximately 212 mm upstream of the 
leading edge of the shock generator. The shock generator itself had a blunt leading edge (not visible 
in the figure). The shock generator extended from approximately x = 173 mm to x = 331 mm, 
and spanned the tunnel. Each of the fine grid zones contained 81 points in the spanwise direction, 
with clustering along the tunnel side wall. The grids were designed to yield minimum spacings 
at walls with y + ~ 1 (approximately 0.005 mm). Near the edge of the incoming boundary layer, 
wall-normal spacing was Ay « 1 mm. 

The flow conditions were: M = 2.25, T ^ = 145.590 K, and Re = 5684 per mm. This 
Reynolds number was based on an assumed freestream density of = 0.1045 kg/m 3 , speed of 
sound of aoo = 241.85 m/s, and viscosity of = 1.0004 x 10 -5 kg/ (ms). Boundary conditions 
were as follows. All tunnel walls (top, bottom, and side) and the shock generator used no- slip 
adiabatic wall conditions. Symmetry was enforced along the centerplane. At the upstream boundary 
located at x = —621 mm, total conditions were specified as: Pt/Poo = 11.56314 and T t /I ^ = 
2.0125. At the downstream boundary located at x = 1000 mm, all quantities were extrapolated 
from the interior of the domain. 

The iterative convergence history of the density residual on the fine grid is shown in Fig. 3. It 
converged a little more than 3 orders of magnitude. The cause of the lack of global convergence 
beyond that point is not known. Sometimes the use of a reconstruction limiter causes convergence 
flattening like this, because of “flip-flopping” (sequentially turning on then off) of the limiter func- 
tion in high gradient areas, e.g., near shock waves. It is also possible that the loss of convergence 
was in part caused by a very large region of separated flow located above and behind the shock 
generator. This separated region was far removed and downstream from the shock boundary layer 
interaction area of interest near the lower wall below the shock generator. The flow field in the area 
of interest did not change noticeably over the course of the final several thousand multigrid cycles. 

A grid study was conducted for this case by running a fine, a medium, and a coarse level. Each 
level was constructed by removing every other grid point from the next finer grid. Results for u- 
velocity profiles and — u'v' profiles at four stations near the interaction region are shown in Figs. 4 


3 



and 5. Overall, there was only a relatively small effect of the grid on ^-velocity, but an exception 
to this was the region near the forward part of the separation zone (x = 300 - 310 mm). In this 
area there were larger changes in ^-velocity due to grid refinement. The overall grid influence on 
—u'v' was even more significant. In this case, the coarse grid was clearly too coarse to adequately 
represent this quantity. The results on the medium and fine grids were closer to each other, again with 
the largest difference near x = 310 mm. The peak —u'v' changed by approximately 6% between the 
medium and fine levels at this location. While this is mostly a qualitative grid sensitivity analysis, 
it gives a feel for the discretization error influence on particular quantities of interest when using 
these grid sizes. It appears that - on the whole - the fine grid may be reasonable for obtaining 
engineering accuracy, but it is probably not fine enough to adequately resolve the region near the 
start of the separation bubble in the interaction region. It would be instructive to learn the influence 
of an extra-fine grid on this sensitive region. 

Centerplane contours of CFD results on the fine grid are compared to experimental PIV data in 
Figs. 6 (a) - (c). The length of the separated region was approximately the same in both the CFD 
and the experiment, although the bubble height was predicted slightly smaller by the CFD. The 
magnitudes of the peak levels of — u'v' were predicted too high by the CFD. However, overall there 
was fairly good correspondence between CFD and experiment in all quantities shown. Profiles at 
specific ^-locations are given in Figs. 7,8, and 9. The ^-velocity profiles were predicted with higher 
magnitude than experiment downstream of x = 300 mm. The u-velocities overall were predicted 
reasonably well, except near x = 320, where the CFD missed the peak near y = 9 mm. This 
difference can also be seen in Fig. 6 (b), where the red contours of high u-velocity in the experiment 
can be seen to extend slightly higher and further downstream than in the CFD solution. 

The three-dimensionality inherent in this case can be seen in the near-bottom wall streamlines, 
shown in Fig. 10. The CFD separated region showed similar features as experiment in the plane 
y = 1 mm above the bottom wall. The CFD bubble in that plane extended from an upstream extent 
near x = 312 mm, to a downstream extent near x = 328 mm, and lateral extent near z & 45 mm. In 
both CFD and experiment there was a curved reattachment line with maximum downstream extent 
at the centerplane. 

4.2 CCAS Cases 

An overall view of the UM fine grid (for 6 = 7.75°) is shown in Fig. 1 1 . All distance units are in mm. 
The grids for the other two CCAS shock angles were similar, so they will not be described in detail. 
However, Fig. 12 shows a comparison of the three shock generator shapes near their leading edges. 
Note that when the grids were created, a very small bluntness of a total height 0.06 mm was specified 
at the leading edges. The grid contained 16 zones, most of which were connected in a point-matched 
(one-to-one) fashion. There were two patched interfaces located at approximately x = — 92 mm and 
at x = 221 mm. The 6 = 7.75° shock generator extended from approximately x = —47 mm to 
x = 121 mm, and spanned only the center region of the tunnel (between 2 = —15.875 to 15.875 mm 
when reflected across the center symmetry plane at z = 0). The shock generator was held in place 
from above by a narrow body located between x = 0 and approximately x = 95 mm. The tunnel 
side wall was located at z = 28.575 mm. The fine grid zones contained 141 points in the span wise 
direction in the region of interest below the shock generator, with clustering along the tunnel side 
wall. The grids were designed to yield minimum spacings at walls with y + ~ 1 (approximately 
0.005 mm). Near the edge of the incoming boundary layer, wall-normal spacing was Ay « 1 mm. 

The flow conditions were: M = 2.75, = 116.617 K, and Re = 8831 per mm. This 

Reynolds number was based on an assumed freestream density of = 0.120 kg/m 3 , speed of 
sound of aoo = 216.4 m/s, and viscosity of = 8.0863 x 10 -6 kg/ (ms). Boundary conditions 
were as follows. All tunnel walls (top, bottom, and side) and the shock generator used no- slip 
adiabatic wall conditions. Symmetry was enforced along the centerplane. At the upstream boundary 
located at x = —949 mm, total conditions were specified as: Pt/Poo = 25.14207 and T t /T ^ = 
2.51251. At the downstream boundary located near x = 600 mm, all quantities were extrapolated 
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from the interior of the domain. 

The global iterative convergence behaviors of the UM cases were worse than the UFAST case. 
The reason for lack of overall convergence over the entire field is not known. However, like the 
UFAST case, the UM cases were run long enough until the flow in the SBLI area of interest near the 
bottom wall did not change appreciably. 

A grid convergence study was performed using two successive grids for the 6 = 7.75° case. The 
^-velocity profiles and —u'v' profiles at four locations along the centerplane in the interaction region 
are shown in Figs. 13 and 14. Overall, the effect of grid (between 0.98 million cells and 7.84 million 
cells) on ^-velocity was relatively small, with the largest influence occurring near the beginning of 
the interaction region, at x = 30 mm. On the other hand, the effect on — u'v' was more significant, 
particularly at the x = 35 mm station. Here, the peak level nearest the wall varied by over 13%, 
and a second peak near y = 5 mm was missed by the coarser grid (this double peak in —u'v' 
occurred underneath the beginning of the reflected shock; there was a similar double peak evident 
in the UFAST experimental data shown in Fig. 9 at x = 320 mm). At other locations, the coarser 
grid underpredicted the peak by as much as about 10%. Thus, the fine grid with 7.84 million cells 
may be reasonable for obtaining engineering accuracy in ^-velocity, but it is not clear how accurate 
it is for turbulence quantities. A finer grid would be required to answer this question. Although not 
readily discernible in Fig. 13, it should be noted that there was no separated flow predicted along the 
centerplane on the medium grid, but a small, very thin region of reverse flow was predicted along 
the centerplane on the fine grid between x = 29 and 36.5 mm. 

Fig. 15 shows the effect of using FDS (as opposed to FVS) on ^-velocities for this case, as well 
as the effect of full Navier-Stokes (as opposed to thin-layer, with viscous derivatives in all coordinate 
directions but no cross-derivative terms). Only the coarser grid with 0.98 million cells was used for 
the results in this figure. FDS had a tendency to yield higher velocity at a given distance from the 
wall than FVS, particularly at x = 30 mm. There was very little difference evident between thin- 
layer and full Navier-Stokes results on the centerplane in the interaction region. Focusing on the 
sensitive x = 30 mm location, the effect of grid density on both FVS and FDS is shown in Fig. 16. 
For both schemes, refining the grid caused higher velocity levels (“fuller” profiles) at this location. 
Although not shown, there was less grid influence at other stations. 

Bottom wall surface-restricted streamlines are shown in Fig. 17. All three methods (FVS thin- 
layer, FDS thin-layer, and FDS full) produced very similar results, including the size and shape 
of the large sidewall interaction region. The main noticeable difference was near the centerplane, 
where the full Navier-Stokes result yielded a very small region of reverse flow, and the others did 
not. 

Contours of CFD results on the fine grid are compared to preliminary experimental PIV data in 
Figs. 18(a) - (b). The PIV grid points (onto which the CFD data were interpolated for the workshop) 
were not close enough to the wall to include any of the very small, thin separated region predicted 
by CFD on the fine grid in this case. Note that there were large inconsistencies in some of the ex- 
perimental data given at the workshop, at the intersections of the span wise planes with the center 
stream wise plane. These inconsistencies were as large as 25 - 35%, and were also noted by DeBonis 
et al. [3], who chose to use only the centerplane measurements for comparisons. Also note that at 
the workshop, some concerns were raised regarding the preliminary experimental turbulence data 
released to the workshop participants. Therefore, no turbulence comparisons will be made here. 
Centerplane CFD comparisons of ^-velocity and u-velocity are shown in Figs. 19 and 20. Predic- 
tions were fairly good for ^-velocity, but were less accurate for u-velocity, although the same trends 
as experiment were exhibited. 

Although not shown, the UM experimental data also included many derivative quantities such as 
du/dx , du/dy , du/dz , etc. However, it it important to note that the PIV grid was used to find these 
derivatives in the experiment. Very near the surface, certain derivatives (such as du/dy) had very 
large gradients - on the order of hundreds of thousands - and they changed rapidly in y. The PIV 
grid density in the ^/-direction (approximately 0.5 mm) was not fine enough to accurately determine 
this particular derivative in the high gradient regions. 
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The blind UM cases 0 = 10° and 9 = 12° were also computed on their respective fine grids. 
Because no experimental data were available for these cases at the time of the workshop, detailed 
results are not given here. However, general CFD-only comparisons between the three UM cases are 
shown in Figs. 21 - 23. The main influence of increasing shock angle was a slightly larger region 
of lower ^-velocity flow in the interaction region on the centerplane. As seen in Fig. 23, all cases 
predicted regions of separated flow near the centerplane, but these regions were very thin so they 
did not show up when the CFD data was interpolated onto the PIV centerplane grids (whose closest 
approach to the wall was around 0.25 mm). The centerplane region was significantly affected by 
very strong sidewall interactions, which were larger for the higher shock angles. 

5 Conclusions 

Numerical simulations of four SBLI cases from the AIAA Supersonic Shock Boundary Layer In- 
teraction Workshop were presented, using the compressible Reynolds-averaged Navier-Stokes CFD 
code CFL3D. The one-equation Spalart-Allmaras turbulence model was employed for all cases. 

Grid convergence studies were conducted for two of the cases. Also, the effect of two different 
flux schemes and the effect of thin-layer vs. full Navier-Stokes was examined. The fine grids (be- 
tween 7 and 9 million cells) were generally dense enough to yield reasonably grid-converged mean 
flow results over much of the interaction region, although the forward part of the interaction region 
on the centerplane was not grid-converged. Turbulence quantities were more grid- sensitive than the 
mean flow overall. Thus, finer grid levels than the ones used for the workshop would be necessary 
to more confidently assess the effects of discretization errors. Flux difference splitting was shown 
to yield somewhat higher velocity levels (“fuller” profiles) than flux vector splitting. The difference 
between thin-layer and full Navier-Stokes on the centerplane in the interaction region was very small 
in the case examined. 

For the UFAST case, CFD predicted centerplane results in the interaction region reasonably well 
compared to experiment, with the exception that the velocity was predicted somewhat too full and 
the peak turbulence levels tended to be overpredicted in magnitude in the region where separated 
flow occurred. For the UM 0 = 7.75° case, ^-velocity comparisons with experiment were overall 
fairly good. The computed u-velocity components exhibited the same trends as experiment, but 
were offset. Comparisons between CFD and experimental — u'v' were not made for the UM case 
because of questions raised during the workshop regarding the latter preliminary data. Strong three- 
dimensional behaviors in and near the interaction regions were documented for all four SBLI cases. 
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Figure 1. Schematic diagram of a SBLI. 
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(c) near x = 190 (centerline symmetry plane on left, side wall on right) 


Figure 2. Views of UFAST fine grid. 
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cycles 


Figure 3. Convergence residual history for UFAST fine grid. 
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Figure 4. Grid effect on ^-velocity profiles on the centerplane in the interaction region of the UFAST 
case. 
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Figure 5. Grid effect on — u'v' profiles on the centerplane in the interaction region of the UFAST 
case. 
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(a) u- velocity, m/s 


(b) ^-velocity, m/s 



(c) — u'v', m 2 /s 2 


Figure 6. Centerplane contours for UFAST case (CFD results from fine grid have been interpolated 
to PIV grid points). 
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Figure 7. Centerplane ^-velocity profiles for UFAST case (CFD results from fine grid have been 
interpolated to PIV grid points). 
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Figure 8. Centerplane v - -velocity profiles for UFAST case (CFD results from fine grid have been 
interpolated to PIV grid points). 
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Figure 9. Centerplane —u'v' profiles for UFAST case (CFD results from fine grid have been inter- 
polated to PIV grid points). 
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(a) CFD bottom wall surface-restricted streamlines, with centerline separation 
extent denoted by arrows 



(b) CFD near y = 1 mm, with centerline separation extent denoted by arrows 



(c) IUSTI experiment near y = 1 mm, from Benek [2] 

Figure 10. Streamlines at and near bottom wall in the interaction region for UFAST case (flow is 
from left to right). 
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(b) along symmetry plane near shock generator 



Figure 11. Views of UM 9 = 7.75° fine grid. 
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Figure 12. Comparison of shock generator shapes for the three different UM cases. 



Figure 13. Grid effect on ^-velocity profiles on the centerplane in the interaction region of the CCAS 
(9 = 7.75° case. 
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Figure 14. Grid effect on —u'v' profiles on the centerplane in the interaction region of the CCAS 
9 = 7.75° case. 
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Figure 15. Effect of FDS and full Navier-Stokes on iz-velocity centerplane profiles in the interaction 
region of the CCAS 6 = 7.75° case, 0.98 million-cell grid. 
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Figure 16. Grid effect on ^-velocity profiles on the centerplane at x = 30 mm for FVS and FDS, 
CCAS <9 = 7.75° case. 
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(a) FVS, thin-layer Navier-Stokes 



(b) FDS, thin-layer Navier-Stokes 



(c) FDS, full Navier-Stokes 


Figure 17. CFD bottom wall surface-restricted streamlines in the interaction region for the CCAS 
6 = 7.75° case, 0.98 million-cell grid (flow is from left to right). 


23 





(a) u- velocity, m/s 



(b) ^-velocity, m/s 


Figure 18. Contours for CCAS 6 = 7.75° case (CFD results from fine grid have been interpolated 
to PIV grid points). 
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Figure 19. Centerplane (/-velocity profiles for CCAS 9 = 7.75° case (CFD results from fine grid 
have been interpolated to PIV grid points). 
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Figure 20. Centerplane ^-velocity profiles for CCAS 6 = 7.75° case (CFD results from fine grid 
have been interpolated to PIV grid points). 
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Figure 21. Comparison of centerplane CFD fine-grid ^-velocity contours for all three CCAS cases 
(CFD results from fine grid have been interpolated to PIV grid points). 
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Figure 22. Comparison of centerplane CFD fine-grid v - velocity contours for all three CCAS cases 
(CFD results from fine grid have been interpolated to PIV grid points). 
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(a) e = 7.75° 



(b) 6 = 10.0° 



(c) 6 = 12.0° 

Figure 23. CFD bottom wall surface-restricted streamlines in the interaction region for all three 
CCAS cases on fine grids (flow is from left to right). 
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